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ABSTRACT 

Gas  flows  in  microfluidic  devices  suffer  from  non-equilibrium  effects.  To  capture  various  non-equilibrium 
phenomena  in  the  slip  and  early  transition  regime,  the  hydrodynamic  system  is  extended  from  the  Navier- 
Stokes-F ourier  equations  to  the  regularised  26  moment  equations  via  Grad’s  moment  method  based  on 
kinetic  theory.  It  is  shown  in  this  lecture  that  the  extended  hydrodynamic  equations  can  be  effectively 
applied  to  micro-electro-mechanical-systems  analytically  and  numerically.  Well  known  non-equilibrium 
phenomenon ,  such  as  Knudsen  layers  and  the  Knudsen  minimum ,  are  used  to  validate  the  extended 
hydrodynamic  model.  It  is  demonstrated  that  the  moment  method  can  be  used  to  predict  gas  flows  in 
micro-electro-mechanical-systems. 


1.0  INTRODUCTION 

Due  to  the  rapid  development  in  fabrication  technology  for  constructing  micro-electro-mechanical- 
systems  (MEMS),  fluid  flow  at  the  micro-  and  nano-scale  has  received  considerable  attention  [1].  A  basic 
understanding  of  the  nature  of  flow  and  heat  transfer  in  these  devices  is  considered  essential  for  efficient 
design  and  control  of  MEMS.  Engineering  applications  for  gas  microflows  include  miniaturised  heat 
exchangers  for  cooling  integrated  circuits,  hand-held  gas  chromatography  systems  for  the  detection  of 
trace  concentrations  of  air-borne  pollutants,  miniature  heat  pumps  and  novel  high-throughput  gas  flow 
cytometers,  etc..  Gas  flows  in  micro-scale  devices  suffer  from  non-equilibrium  effects  when  the  gas 
molecular  mean  free  path  is  the  same  order  as  the  characteristic  length  of  the  device.  The  degree  of  non¬ 
equilibrium  of  a  gas  is  generally  expressed  through  the  Knudsen  number  ( Kn=A/L )  which  is  the  ratio  of 
the  molecular  mean  free  path,  A,  to  a  typical  dimension  of  the  flow  field,  L.  The  non-equilibrium  gas  flow, 
or  rarefied  gas  dynamics  has  been  explored  extensively  for  more  than  a  century  in  association  with  high¬ 
speed  high-altitude  flow  applications,  such  as  space  re-entry  vehicles,  and  flows  under  ultra-low  pressure 
(vacuum)  conditions,  where  A  has  a  large  value.  Gas  flow  in  micro-scale  devices  can  suffer  from 
rarefaction  effects  because  the  characteristic  length  of  the  device,  L,  is  so  small  that  it  is  comparable  to  the 
mean  free  path  of  the  gas,  even  under  atmospheric  conditions.  Beskok  [2]  compiled  a  typical  range  of  gas 
flow  regimes  in  terms  of  micro  device  dimension,  as  show  in  Figure  1.  With  the  advent  of  MEMS,  there 
has  been  a  renewed  impetus  in  the  development  of  new  and  efficient  approaches  for  modelling  low-speed 
slip-  and  transitional-flows  that  can  capture  non-intuitive  phenomena,  such  as  Knudsen  layers  and  the 
Knudsen  minimum  [3],  and  provide  an  accurate  description  of  a  gas  that  is  not  too  far  from 
thermodynamic  equilibrium  [4,  5]. 
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Figure  1:  Characteristic  length  scales  of  typical  microfluidic  components  and  the  corresponding 
Knudsen  number  at  standard  atmospheric  conditions.  Adapted  from  Beskok  [2]. 


The  key  flow  features  for  MEMS  are  that  flow  is  laminar,  typically  Reynolds  number  Re  <  1  and  Mach 
number  Ma  <  1  and  geometry  can  be  complex.  The  different  rarefaction  regimes  can  be  summarised 
according  to  the  value  of  the  Knudsen  number  [6]:  (i)  no  slip  (Kn  <  10'3);  (ii)  slip  (10'3<  Kn  <10_1);  (iii) 
transition  (10'1  <  Kn  <  10);  and  (iv)  free  molecular  flow  Kn  >  10,  respectively.  Most  MEMS  operate  at  the 
slip  and  early  transition  regime  (Kn  <  1 )  as  shown  in  figure  1 .  Different  approaches  have  been  employed 
by  various  researchers  to  capture  and  describe  the  non-equilibrium  phenomena  that  arise  due  to  an 
insufficient  number  of  molecular  collisions  occurring  under  rarefied  conditions.  These  approaches  can  be 
broadly  divided  into  two  specific  categories:  the  microscopic  and  the  macroscopic.  In  this  lecture,  the 
emphasis  will  be  on  the  development  of  the  macroscopic  methods  and  their  application  towards 
engineering  configurations.  Microscopically,  the  Boltzmann  equation  [7]  provides  an  accurate  description 
of  a  dilute  gas  at  all  degrees  of  rarefaction  and  describes  its  state  through  a  molecular  distribution  function 
that  treats  the  gas  as  a  large  number  of  interacting  molecules,  colliding  and  rebounding  according  to 
prescribed  laws.  However,  solutions  of  the  Boltzmann  equation,  either  directly  [8,  9]  or  through  the  direct 
simulation  Monte  Carlo  (DSMC)  method  [10],  entails  significant  mathematical  complexity  and  can  be 
computationally  expensive,  particularly  for  low-speed,  low  Knudsen  number  flows  in  the  slip  and 
transition  regime. 

Due  to  the  difficulties  associated  with  solving  the  Boltzmann  equation,  there  is  significant  effort  being 
made  to  construct  alternative  solution  strategies  that  can  provide  an  accurate  description  of  a  gas  with 
Knudsen  numbers  that  extend  into  the  transition  regime.  For  designing  components  in  MEMS,  it  is 
desirable  that  any  new  developments  have:  (i)  computational  efficiency  comparable  to  conventional 
hydrodynamic  formulations;  (ii)  the  ability  to  handle  real  geometries,  and  (iii)  under  appropriate 
conditions,  will  recover  the  Navier-Stokes-Fourier  solution.  The  two  main  approaches  used  to  derive  these 
extended  hydrodynamic  (EHD)  equations  are  the  Chapman-Enskog  (CE)  expansion  [11]  and  the  method 
of  moments  developed  by  Grad  [12], 

The  Chapman-Enskog  approach  expands  the  molecular  distribution  function  in  powers  of  Kn  to  construct 
the  constitutive  relationships  using  the  conventional  hydrodynamic  variables  and  their  derivatives.  The 
zeroth-order  expansion  yields  the  Euler  equations  and  the  first-order  results  in  the  Navier-Stokes-Fourier 
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(NSF)  equations.  In  association  with  velocity-slip  and  temperature-jump  wall  boundary  conditions,  the 
NSF  equations  can  predict  the  velocity  profile  in  a  micro-channel  fairly  accurately  in  the  slip  regime  for 
isothermal  flows.  Higher  order  expansions  yield  the  Burnett  equations  (second-order),  Super-Burnett 
equations  (third-order),  and  so  on,  which  can  be  unstable  in  time  dependent  problems.  In  Grad’s  approach, 
a  set  of  governing  partial  differential  equations  representing  the  moments  of  the  molecular  distribution 
function  are  derived  from  the  Boltzmann  equation.  However,  moments  of  higher  order  always  appear  in 
each  moment  equation  and  the  set  of  moment  equations  is  not  closed.  To  avoid  the  necessity  of  dealing 
with  an  infinite  number  of  moment  equations,  a  closure  procedure  is  required  that  relates  the  higher-order 
moments  to  those  of  lower  order.  In  the  seminal  work  of  Grad  [11],  the  set  of  moment  equations  was 
closed  at  the  second-moment  level,  which  involves  13  moments:  density,  momentum,  energy,  heat  flux, 
and  pressure  deviator.  It  is  interesting  to  point  out  that  the  constitutive  relationships  established  by  the 
Chapman-Enskog  method,  at  any  order,  can  be  regarded  as  a  first  moment  closure  method  with  5 
moments.  To  close  the  set  of  moment  equations  at  the  second-moment  level,  Grad  [11]  expanded  the 
distribution  function  in  Hermite  polynomials  about  the  local  Maxwellian  to  third-order  accuracy  and  set 
the  trace-free  part  of  the  third  moments  to  zero.  As  a  result,  the  original  13  moment  equations  derived  by 
Grad  (G13)  are  hyperbolic  and  lack  any  gradient  transport  mechanism.  They  are  not  suitable  for 
computing  boundary  layers  which  appear  in  MEMS. 

Struchtrup  and  Torrilhon  [13]  regularised  Grad’s  13  moment  equations  (R13)  by  applying  a  CE-like 
expansion  to  the  governing  equations  of  the  moments  higher  than  second  order.  Algebraic  constitutive 
relationships  were  then  established  between  the  higher  moments  and  the  second  and  lower  moments.  The 
R13  equations  can  give  reliable  results  for  Knudsen  numbers  up  to  Kn  <  0.25.  To  extend  the  capability  of 
the  moment  method  to  larger  Knudsen  numbers  in  the  early  transition  regime,  we  extended  the 
regularisation  theory  to  the  26  moment  equations  (R26)  [14].  However,  to  apply  the  regularised  moment 
equations  to  confined  flows,  such  as  those  found  in  microfluidic  channels,  wall  boundary  conditions  are 
required  for  higher  moments  in  addition  to  the  velocity-slip  and  temperature-jump  conditions.  We 
pioneered  the  treatment  of  wall  boundary  conditions  for  macroscopic  variables  from  kinetic  boundary 
conditions  [15].  The  treatment  of  solid  boundaries  was  further  improved  by  Torrilhon  and  Struchtrup  [16]. 
This  opens  the  door  to  solve  the  moment  equations  for  flows  in  confined  geometries. 

For  simple  geometries,  such  as  a  straight  channel  or  pipe,  analytical  solutions  can  be  obtained  from  the 
linearised  moment  equations  [17-21].  However,  for  real  engineering  applications  of  complicated 
geometries,  numerical  procedures  have  to  be  employed  to  obtain  solutions.  To  numerically  solve  the 
moment  equations  for  parabolic  and  elliptic  flows  within  a  conventional  computational  fluid  dynamics 
(CFD)  approach,  such  as  the  finite  volume  (FV)  method,  is  quite  challenging  because  the  gradient 
transport  mechanisms  are  not  explicitly  expressed  in  the  momentum  and  energy  equations  of  the  moment 
system.  The  inadequacy  of  the  standard  FV  method  for  the  governing  equations  without  any  gradient 
transport  terms  is  well  recognised  [22],  so  that  methods  for  dealing  with  hyperbolic  problems  are  required, 
such  as  Riemann  solvers  [23],  TVD  schemes  [24]  and  ENO  schemes  [25].  These  schemes  are  complex 
and  computationally  expensive  for  multidimensional  confined  flow,  particularly  at  low  speed.  Gu  and 
Emerson  [14,  15]  proposed  a  primitive  variable  transformation  approach  to  split  the  higher  order  moments 
into  their  gradient  and  non-gradient  transport  components.  Recently,  a  new  strategy  is  proposed  to  solve 
the  regularised  moment  equations  on  a  collocated  grid  for  low  speed  flows  in  confined  geometries  [26],  In 
this  lecture,  the  approach  to  build  up  macroscopic  governing  equations  by  the  moment  method  will  be 
introduced.  The  procedure  to  solve  the  moment  equations  for  low  speed  flows  in  MEMS  will  be  described 
and  validation  examples  will  be  given. 


2.0  MOMENT  METHOD 

Kinetic  theory  accounts  for  a  molecule’s  movement  and  interaction  through  a  molecular  phase-density 
distribution  function,  F( £  x,  t),  which  satisfies  the  Boltzmann  integro-differential  equation  [6],  where  x 


RTO-EN-AVT-194 


11  -3 


Application  of  the  Moment  Method  in  the 

Slip  and  Transition  Regime  for  Microfluidic  Flows 


ORGANIZATION 


and  4  arc  the  position  and  velocity  vectors,  respectively,  of  a  molecule  at  time  t,  and  F(4  x,  t)dxd%  gives 
the  number  of  molecules  whose  velocities  lie  within  d%  in  a  volume  element  dx.  For  convenience,  a  mass 
distribution  function  is  used  in  the  present  study  and  is  defined  by 

f(£,x,t)  =  mF(%,x,t),  (1) 


where  m  is  the  mass  of  a  molecule.  The  Boltzmann  equation,  expressed  by  [6]: 


dJL 

dt 


,  df  df 

+  —  +  a  ■ - : 

d$ 


Q{fj) 


(2) 


is  the  central  equation  in  kinetic  theory,  the  properties  of  which  can  be  used  to  guide  the  development  of 
kinetic  and  macroscopic  models  for  rarefied  gas  flow.  Here  t  and  x,  are  temporal  and  spatial  coordinates, 
respectively,  and  any  suffix  i,j,  k  represents  the  usual  summation  convention.  The  external  acceleration  is 
denoted  by  a,-.  On  the  right  hand  side  of  equation  (2),  Q(f,  f)  is  the  collision  integral.  The  most  interesting 
feature  of  the  Boltzmann  equation  is  that  the  macroscopic  irreversible  process  of  entropy  is  imbedded  in  it, 
known  as  the  //-theorem  [6], 

Kinetic  theory  associates  hydrodynamic  quantities  with  the  moments  of  the  distribution  function.  Once  the 
distribution  function,  f  is  known,  its  moments  with  respect  to  £,  can  be  determined.  For  example,  the 
density,  p ,  and  the  momentum,  pa,,  can  be  obtained  from  [27]: 

P  =  \fd4  and  PUi=\4ifd4.  (3) 


An  intrinsic  or  peculiar  velocity  is  introduced  as 


ci  =  £i~ui>  (4) 

so  that  the  moments  with  respect  to  «,  can  be  conveniently  calculated  as 

Phh . A  =lchch . £tJdl  (5) 

The  pressure  tensor  is  a  second  moment  of  the  distribution  function  and  can  be  separated  as  follows: 

Pij  =  |  Cfijf  d£  =  P8ij  +  P<ij>  -  pSy  +  Gy  (6) 


where  Sy  is  the  Kronecker  delta  function,  p  =  pkkJ 3  is  the  pressure,  and  Gy  =  p<y>  is  the  deviatoric  stress 
tensor.  The  angular  brackets  denote  the  traceless  part  of  a  symmetric  tensor.  Furthermore,  the  thermal 
energy  density  and  the  heat  flux  vector  are  given,  respectively,  by 

ps  =  Uc2fd4  =  ^p-T  =  ^pRT,  (7) 

2  2  m  2 


di=\\c2  cif  d4  ■  (8) 

where  c2  =  ckck-  The  temperature,  T,  is  related  to  the  pressure  and  density  by  the  ideal  gas  law, 
p=p{klm)T=pRT,  where  k  is  Boltzmann's  constant  and  R  is  the  gas  constant. 
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Equation  (2)  tells  us  that  there  is  only  one  equilibrium  state  and  its  distribution  function  is  expressed  by 
the  Maxwellian,  fM,  parameterised  by  traditional  hydrodynamic  variables,  p,  T,  and  u,  as  [6]: 


fit  ~ 


P 


{In  RT) 


-exp 


2  RT 


(9) 


When  the  gas  is  away  from  the  equilibrium  state,  an  infinite  number  of  independent  moments  can  be 
defined  from  equation  (5).  In  other  words,  to  recover  the  distribution  function,  f  at  the  state  of  non¬ 
equilibrium,  an  infinite  number  of  moments  are  required.  In  the  moment  method,  a  limited  number  of 
moments  are  employed  to  approximate  the  non-equilibrium  state. 

The  five  collision  invariants,  y/=  (1,  c2),  which  satisfy  [6]: 

WQ{f,fY$  =  0,  (10) 


lead  to  the  conservation  laws  for  mass,  momentum  and  energy,  respectively,  as  [27] : 

dt  dx[ 


(11) 


dpi!:  dp  up  dati 

- -  + - -  +  — - 

dt  8xj  dxj 


dp 
dx i 


+  Pa i 


dpT  dpiijT  2  dql 
dt  dxt  3  R  dxl 


2  f  diL  du  , 

—  P — L  +  <Ja — ~ 

3  R  ^  dXj  1  dxi  J 


(12) 


(13) 


To  close  the  above  set  of  equations,  a  CE  expansion  of  the  distribution  function  in  powers  of  Kn  is  applied 
to  the  Boltzmann  equation.  The  traditional  hydrodynamic  variables,  p,  T,  and  u,  and  their  derivatives  are 
used  to  estimate  the  values  of  Gy  and  q,.  Alternatively,  the  governing  equations  for  Gy,  q,  and  higher 
moments  can  be  derived  from  the  Boltzmann  equation.  However,  there  appear  even  higher  order  moments 
in  the  relevant  moment  equations.  There  are  various  ways  to  close  the  moment  equation  set  [28]. 

In  the  moment  method  proposed  by  Grad  [11],  the  distribution  function,  /  is  expanded  in  Hermite 
polynomials  as: 


/  =  I™  fu  X  —,a(AH(A]  -  fit 


n=0  n\ 


a(0)H(0)  +  +  ^af  Hf  +  + .. 


(14) 


where  H\  ’  is  the  Hermite  function  and  a\  ’  are  the  coefficients.  All  of  the  polynomial  coefficients  are 
linear  combinations  of  the  moments  of  f.  An  infinite  set  of  Hermite  coefficients  is  equivalent  to  the 
molecular  distribution  function  and  no  kinetic  information  is  lost  in  such  a  system.  In  practice,  the 
molecular  distribution  function  has  to  be  truncated  and  the  specific  problem  to  be  addressed  will  determine 
the  order  of  the  truncation.  Usually,  the  truncated  distribution  function  is  called  Grad’s  distribution 
function  and  is  denoted  as  /c. 

It  is  expected  that  as  Kn  increases,  the  order  of  the  truncated  Hermite  polynomials  should  increase  to 
provide  an  accurate  description  of  the  flow  conditions.  By  choosing  a  sufficient  number  of  Hermite 
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coefficients,  a  general  solution  of  the  Boltzmann  equation  can  be  approximated.  The  governing  equations 
of  the  moments  involved  in  the  Hermite  coefficients  can  be  derived  from  the  Boltzmann  equation. 
However,  the  set  of  equations  are  not  closed,  since  the  fluxes  and  the  collision  terms  are  unknown.  The 
truncated  distribution  function  is  used  to  calculate  the  collision  terms  and  higher  moments  in  the  fluxes  as 
functions  of  the  chosen  moments  to  close  the  equation  set.  All  the  moments  included  in  the  truncated 
distribution  function  construct  the  Grad  moment  manifold  (GMM).  These  moments  relax  to  the 
equilibrium  state  at  a  rate  governed  by  their  corresponding  governing  equations.  In  Grad’s  original 
approach,  the  remaining  higher  moments  outside  the  GMM  are  calculated  as  a  linear  combination  of  the 
moments  in  the  GMM  without  any  dynamic  corrections.  In  the  regularised  moment  method,  a  CE 
expansion  is  applied  to  the  higher  order  moment  equations,  so  that  the  remaining  higher  moments  outside 
the  GMM  from  the  truncated  distribution  function  approach  the  GMM  (not  the  equilibrium  state)  at  a  fast 
finite  rate,  then  relax  to  the  equilibrium  state  along  with  the  GMM  [28].  The  different  relaxation  process 
models  used  in  the  construction  of  macroscopic  equations  are  schematically  shown  in  Figure  2.  Circles 
with  broken  lines  represent  a  non-equilibrium  state  of  an  infinite  number  of  independent  moments.  Circles 
with  solid  line  represent  an  equilibrium  or  quasi-equilibrium  sate  which  can  be  represented  by  a  finite 
number  of  independent  moments.  It  clearly  illustrates  the  essential  role  of  the  CE  expansion  in  the 
modelling  process.  An  infinite  number  of  independent  moments  can  be  absorbed  into  a  system  with  a 
finite  number  of  moments  through  a  CE  expansion. 


non-equilibrium  equilibrium  quasi-equilibrium  equilibrium 

t'  \  NSF equations 

/  , - m - J/m) 

\  /  CE  expansion 

(a)  CE  expansion  on  the  Boltzmann  equation  (b)  Grad’s  moment  method 


non-equilibrium  quasi-equilibrium  equilibrium 


(c)  Regularised  moment  method 

Figure  2:  Schematic  modelling  of  relaxation  process  in  the  moment  method. 


2.1  Continuum  Hydrodynamic  Model  -  Navier-Stokes-Fourier  Equations 

The  traditional  hydrodynamic  quantities  of  density,  p,  velocity,  uh  and  temperature,  T,  correspond  to  the 
first  five  lowest-order  moments  of  the  molecular  distribution  function.  The  governing  equations  of  these 
hydrodynamic  quantities  for  a  dilute  gas  can  be  obtained  from  the  Boltzmann  equation  and  represent  mass, 
momentum,  and  energy  conservation  laws,  as  expressed  in  equations  (1 1)-(13).  The  classical  way  to  close 
this  set  of  equations  is  through  a  CE  expansion  of  the  molecular  distribution  function  in  terms  of  Kn 
around  the  Maxwellian,  which  is  first  order  in  Hermite  polynomials.  The  zeroth-order  CE  expansion 
yields  the  Euler  equations  and  the  first-order  approximation  of  cry  and  q„  for  Maxwell  molecules,  gives 
[10,27]: 


aNSF=_2 

ox}> 


.  nsf  15  dT 
and  qt  =  -—Rp—, 

4  OX; 


(15) 


in  which  p  is  the  viscosity.  Equation  (15)  expresses  an  import  transport  mechanism  for  u,  and  T:  the 
gradient  transport  mechanism  in  continuum  mechanics.  If  we  let 
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o,,  =  cr, 


NSF 


and  qt=qt 


NSF 


(16) 


and  inserting  equations  (15)  and  (16)  into  equations  (12)  and  (13)  results  in  the  traditional  hydrodynamic 
equations,  i.e.  the  Navier-Stokes-Fourier  equations.  The  second  order  CE  expansion  adds  the  higher 
derivatives  of  the  hydrodynamic  variables  to  cfy  and  q,  and  results  in  the  Burnett  equations. 


2.2  Second  Moment  Closure  -  Regularised  13  Moment  Equations 

As  the  value  of  Kn  increases,  more  moments  need  to  be  included  in  the  GMM  to  accurately  describe  any 
non-equilibrium  phenomena.  Grad  [11]  truncated  the  distribution  function  to  the  incomplete  third  order  in 
Flermite  polynomials,  fan,  as: 


fc  13  -  ft 


M 


2  pRT 


c  c  ■  — 

‘  j 


pRT 


cicli 


1  — 


5  RT 


(17) 


Grad  was  one  of  the  pioneers  to  introduce  (Jy  and  q,  as  extended  hydrodynamic  variables  and  derived  a  set 
of  governing  equations  for  them  from  the  Boltzmann  equation.  For  Maxwell  molecules,  the  stress  and  heat 
flux  equations  are  [27]: 


do. 


8t 


ij  ,  8uk°ij  , 


dx. 


dx, 


P 


8Xj> 


4  oq ,  dU  j  > 

■  ~  2 cr. 


5  dxf> 


k< 


dx. 


(18) 


8q  duJqi  \  GRlf 
— —  +  — - —  + - - 

dt  dxj  2  dxj 


3  P 


dx , 


2  dx,. 


dx,. 


7  du,  du,  du. 


dx,. 


dx, 


dx, 


1  dA 


P 

du 


dp 

KdxJ 


do 


jk 


dx , 


k  7 


- miik  — - 

6  dx,  dx,. 


(19) 


If  only  the  underlined  terms  in  equations  (18)  and  (19)  are  retained  and  the  rest  of  terms  are  set  to  zero,  the 
NSF  equations  are  recovered.  Flere,  mijk,  Rj  and  A  represent  the  difference  between  the  true  value  of  the 
higher  moments  ( p<jk> ,  p<y>m  and  Pms)  and  their  approximated  value  with  fG13,  respectively  [13],  i.e. 


^hjk  —  P<ijk>  Fijk>\fGi2  ~  P<ijk>’ 

Fj  ~  P<ij>rr  P<ij>rr\fg\3  ~  P<ij>rr  ^RTOy,  (20) 

^  —  Prrss  Prrss\fgi3  Prry,  15  pRT . 


O’]  q  O'  I  q  o'  i  q 

In  Grad’s  original  method,  such  deviations  were  omitted,  so  that  mijk  =Rij  =A  =0.  This  results  in 

the  well  known  G13  equations.  To  close  the  set  of  equations  (1 1)— ( 13),  (18)  and  (19),  Struchtrup  & 
Torrilhon  regularised  the  G13  equations  and  obtained  the  following  closures  [13,  27]: 


mS3=~2  ~ 
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RT- 


do 


<y 


dx, 


+  o 


dRT 


<j 


dp 


<y 


k> 


dx 


k> 


P  dxk> 


4 

+~P 


dU  : 


dx, 


k> 


a<v  Sok>, 
p  dx, 


(21) 
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If  we  let 


r>fil3  24  // 

»  ='T  ~P 


RTd^L  +  2qJRT  “« 


dx 


i> 


dx 


j> 


dp  +  daj>k 


\dxj> 


dx. 


,5f£ 

6  p 


dqk 

ydxk 


+  cr. 


du± 

dx, 


i 


5RT 


r 


k<i 


du j> 

dx. 


+  cr 


du,. 


k<i 


dx 


J> 


2  du. 


4  ®k<i® j>k 
7  p 


(22) 


A^13  =  —12— 


7  srt  q, 
- qk -z - - 


P 


dp_+doJL 
ydxj  dxk 


+  RT 


dqk 


du 


ydxk  U'1dxJJ 


(23) 


mijk=mijk  ,Ry=Ry  and  A  =  A 


R13 


(24) 

and  inserting  equations  (21)-(24)  into  equations  (18)  and  (19),  we  obtain  a  set  of  closed  13  moment 
equations  which  is  denoted  as  the  R13  equations  [13]. 


2.3  Higher  Order  Moment  Closure  -  Regularised  26  Moment  Equations 

One  of  the  challenges  for  macroscopic  models  in  microflows  is  to  capture  the  Knudsen  layer  [4],  The  R13 
equations  improved  Grad’s  original  closure  significantly,  but  they  were  not  able  to  capture  the  Knudsen 
layer  velocity  profile  accurately  [21].  Since  equations  (21)  and  (22)  are  algebraic  approximations  for  mi;t 
and  Rij,  they  can  only  provide  a  mechanism  to  produce  a  boundary  layer  for  the  lower  order  moments,  a,, 
and  q„  but  have  no  mechanism  to  produce  their  own  boundary  layer  near  the  wall.  As  the  Knudsen  layer  is 
a  linear  supeiposition  of  many  exponential  layers  in  planar  [29],  the  R13  equations  only  resolve  one  such 
contribution  [21],  To  increase  the  capability  of  the  moment  method,  the  distribution  function  is  truncated 
to  the  incomplete  fourth  order  in  Hermite  polynomials  as  [14]: 


fa 26  ~  Im  v  + 
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f  „2 
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f  c2 

1RT 


-1 


15(RT) 


3  RT 


+  1 


/  J 


(25) 


and  the  moments  mljk,  R,,  and  A  are  included  in  the  GMM  as  extended  hydrodynamic  variables.  As  a  result, 
a  set  of  26  moment  equations  can  be  derived  from  the  Boltzmann  equation  [14]: 
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(27) 
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SA  dA.lL  oQ, 
8t  8xl  dx. 


3// 


dRT 

dxk 


clz 

f  „  -  X 

dp  |  dajk 

+  RT 

dqk  !  ^  dut  "1 

p 

{ dxj  dxk  J 

[dxk  lJ  dx j) 

2  RT(7,j(7 ij 

3  P 


.(28) 


Here,  y/tjk  and  £2;  are  the  difference  between  the  true  value  of  the  higher  moments  (p<yk> ,  p<ijk>rr,  and 
p,TSsi)  and  their  corresponding  value  approximated  with  fG26,  i.e. 

fiijkl  —  P<ijkl>  P<ijkl> \fG26  ~  P<ijkl>  ’ 

y  ijk  —  Prr<ijk>  Prr ijk >  ~  Prr<ijk>  9RTm \yk ,  (29) 

=  Prrssi  ~  Prrssi\fG26  =  Prrssi  ~  3SRTqt. 


Following  the  same  regularisation  procedure  for  the  13  moment  equations,  the  following  closures  for 
equations  (26)-(28)  can  be  readily  obtained  by  [14]: 
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(  5wpo> +  14^/7/ 

15 
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(32) 


in  which  Zi=2.097,  Z2=0.291,  7i=1.698,  72=l-203  and  73=0.854  are  collision  term  constants  for  Maxwell 
molecules  [14,  30].  Following  the  convention  of  Struchtrup  [27],  the  above  closed  set  of  26  moment 
equations  are  denoted  as  the  R26  equations.  In  the  right  hand  side  of  the  governing  equations  (26)-(28)  for 
the  moments  mijk,  R,,  and  A,  only  the  terms  involved  in  the  R13  closure  are  retained.  Some  higher  order 
terms  in  the  constitutive  relationships  (30)-(32),  which  have  negligible  effects  for  flows  in  MEMS,  are 
neglected.  The  full  list  of  equations  and  their  closures  can  be  found  in  Ref.  [14]. 


4.0  WALL  BOUNDARY  CONDITIONS 

To  apply  any  of  the  foregoing  models  to  flows  in  confined  geometries,  appropriate  wall  boundary 
conditions  are  required  to  determine  a  unique  solution.  One  of  the  difficulties  encountered  in  any 
investigation  of  wall  boundary  conditions  is  due  to  a  limited  understanding  of  the  structure  of  surface 
layers  of  solid  bodies  and  of  the  effective  interaction  potential  of  the  gas  molecules  with  the  wall.  A 
scattering  kernel  represents  a  fundamental  concept  in  gas-surface  interactions,  by  means  of  which  other 
quantities  should  be  defined  [2].  For  the  moment  equations,  the  boundary  conditions  for  the  relevant 
macroscopic  variables  are  required  in  addition  to  the  traditional  velocity  slip  and  temperature  jump 
conditions.  As  we  can  derive  the  macroscopic  governing  equations  from  kinetic  theory,  the  macroscopic 
boundary  conditions  can  also  be  constructed  from  kinetic  theory.  Gu  &  Emerson  [14,  15]  obtained  a  set  of 
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wall  boundary  conditions  for  the  R13  moment  equations  based  on  Maxwell’s  kinetic  wall  boundary  model 
[31]  and  a  fourth-order  approximation  of  the  molecular  distribution  function  in  Hermite  polynomials.  To 
construct  wall  boundary  conditions  for  the  R26  equations,  a  fifth  order  approximation  of  the  molecular 
distribution  function,  /(:,) ,  is  required,  which  can  be  expressed  as  [14]: 


/(  ^  -  fc  26  +  Af 


rijkicicjckci 


+ 


Yijkcffk 


24  p(RT)  12  p(RT) 


r  c2 

9RT 


-1 


CA 


40  p(RT) 


7  (RT) 


2c? 

RT 


+  5 


(33) 


It  should  be  noted  that  the  higher  moments  involved  in  the  underlined  terms  in  equation  (33)  arise  from 
the  regularisation  procedure. 

Maxwell's  kinetic  boundary  condition  [31]  is  one  of  the  simplest  models  and  it  states  that  a  fraction,  a, 
will  be  diffusely  reflected  with  a  Maxwellian  distribution,  ff ,  at  the  temperature  of  the  wall,  Tw,  while  the 
remaining  fraction  (1  -a),  of  gas  molecules  will  undergo  specular  reflection.  In  a  frame  where  the 
coordinates  are  attached  to  the  wall,  with  n,  the  normal  vector  of  the  wall  pointing  towards  the  gas  and  r, 
the  tangential  vector  of  the  wall,  such  that  all  molecules  with  <  0  are  incident  upon  the  wall  and 
molecules  with  >  0  are  emitted  by  the  wall,  Maxwell’s  boundary  condition  can  be  expressed  by  [27]: 


r= 


\  afM  +(!-«)  fi-f'h  ),  >  0, 


/(£«,), 


Zini  <  0. 


(34) 


By  definition,  the  value  of  any  moment  at  the  wall  can  be  obtained  from 

Wo chch  ••• cin  +  (l-«)/K",)]^.  (35) 

From  the  analysis  of  the  special  case  of  a  =  0,  it  was  found  that  only  moments  that  are  odd  functions  of 
cinl  can  be  used  to  construct  wall  boundary  conditions  [11].  Furthermore,  only  those  moments  representing 
fluxes  can  be  used  in  the  boundary  conditions  [16].  From  equation  (35),  the  slip  velocity  parallel  to  the 
wall,  u„  and  temperature  jump  conditions  can  be  obtained  as  [14,  15]: 


a 


u  2-q  jvRT  a m  5mmn  +  2qT  |  9Qr  +  70i//„ 
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l0/T 


2520 paRT 


(36) 


RT-RT=  — 
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2- a  \kRT  qn  RTan 


2  2P  a  4Pa 


-  +  -^~ 
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75R„„+28A 
840 Pa 


24  Pa 


(37) 


where 


p  =  p  +  ^JHL-30R^+7A  (38 

2  840 RT  24  RT 

Here  crnn  (Tjtiifij ,  cy:! T  cTijii i Zj.  q z  q : p,  q n  cj ,h i.  iwnnT  in jqi iHj Tk*  utnnn  iw  jiji qijii Run  Ryu Pj.  T1t  f -> r,, 

y/nn t  ~  WijinPjTk,  and  (j),mnn  =  fikinp/iqii  are  the  tangential  and  normal  components  of  cry,  qt,  mijk,  Ry,  y/yk,  Q 
and  (pyki  relative  to  the  wall,  respectively.  It  should  be  noted  that  the  normal  velocity  at  the  wall,  un  =  0, 
since  there  is  no  gas  flow  through  the  wall.  Equations  (36)  and  (37)  are  similar  to  the  slip  velocity  and 
temperature  jump  conditions  for  the  NSF  equations  [5,  6]  with  the  underlined  terms  on  the  right  hand  side 
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providing  higher-order  moment  corrections.  These  underlined  terms  can  be  related  to  second-  or  higher- 
order  velocity  slip  and  temperature  jump  boundary  conditions.  The  slip  velocity  in  the  higher  moment 
system  is  not  only  proportional  to  the  shear  stress,  cr„„  but  also  to  the  tangential  heat  flux,  q„  and  other 
higher  moments.  With  a  normalised  slip  velocity,  uT  —  ut/\JrT  ,  and  a  normalised  wall  temperature, 
Tw  =TW/T ,  the  rest  of  the  wall  boundary  conditions  are  [14]: 
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5.0  LINEARISATION  OF  MOMENT  EQUATIONS 

One  of  the  main  features  for  gas  flow  in  MEMS  is  low  Reynolds  number,  Re  <  1  and  low  Mach  number 
Ma  <  1.  In  most  cases,  it  is  sufficient  to  apply  a  linearised  moment  equation  set  to  flow  in  the  micro¬ 
system.  For  flows  with  small  deviations  from  an  equilibrium  state  given  by  pa,  Ta,  (or,  p„)  and  uii0,  these 
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deviations  can  be  represented  by  dimensionless  variables  p,  T,  w(,  cr-,  qt,  mijk,  RP  A > 
$ijki>  Vijk  ancl  Q,  in 


P  =  Po{l  +  P )> 

°H=PoRTo°ip 

A  =  p0{RT0fA, 

xi=Lxi, 


t  =  t0(\+t), 

Qi  =  Po  iRRo)31  Qi’ 
<t>ijkl=  Po(RTof  frjkb 
t  =  (L/jRTa)T, 


P=Po{l  +  P )> 

mijk=  Po(RTof2  ™ijk > 
Yijk  =  Po{RTof12  Wyk, 


ui  =ylR^ui, 

Rij  =  Po{RTo)2  Rij’ 
n,=Po(RT0f2hi, 


(47) 


in  which,  L  is  a  characteristic  length  of  the  microsystem.  Inserting  equations  (47)  into  equations  (11)  to 
(13),  (1 8)-(  1 9)  and  (26)-(28)  and  keeping  only  the  terms  that  are  linear  in  the  deviations,  a  set  of  linearised 
R26  moment  equations  (LR26)  are  obtained  as  [19]: 


3£+Al  =  0, 

dt  dxf 

(48) 

dii:  da-j-  dp  La 

_L  + - A- - £_  + - L5 

Sx)  3x)  RZ) 

(49) 

dp  2  dqt  5 

cZ  3  ox)  3  ax)  ’ 

(50) 

5^  |  owi/t  _  p  Z  _  o  0m<(.  4  5g<; 

dt  dxk  \  2  A  >j  dxj>  5  dxj>  ’ 

(51) 

dq,  1  oRy  2  pf_  5  SR  daik  1  5A 

dt  2  dxj  3\  2  A^'  2  dxi  dxk  6  dxi 

(52) 

8mijk  _  3  IJl  dtr<(/ 

dt  dx,  2  V  2  A  iJk  dxk> 

dRij  |  dVijk  _  7  I^L-  28  dq<t 
dt  dxk  6  v  2  A  ij  5  dXj> 


(53) 

(54) 


dA  dQ ,  2  Ur  L  -  adq; 

—  + - L  = —  J - A -8—. 

dt  dx[  3  v  2  2  &,■ 


The  linearised  constitutive  relationships  are: 

t  f2~ T  dira^ 

'yi/  Zx\nL  dxly 


(55) 


(56) 
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Yijk 


27  [2  A  dR 


<y 


lYx\n  L  dxj 


k> 


2  A  dR, ■, 


3  V  n  L  dx,  \  n  L  dx, 


(57) 

(58) 


along  with  the  linearised  boundary  conditions: 
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r  a  V2  m  10  2520 


(59) 


T_  2~«  X  o„n  75Rnn  +  28 A  !  1 

840 


a  \  2  2  4 


h  T 

nnnn  _|_  w  _ ^ 


24  T 


cr,  =  — 


2 -a  \n 
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(60) 
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(66) 

(67) 

(68) 


in  which, 


A  =  - 


Po 


M  \nRT0 


(69) 


is  the  mean  free  path.  For  flows  in  simply  geometries,  such  as  Couette  and  Poiseuille  flows,  the  above 
linearised  moment  equations  can  be  decoupled  into  three  sets  of  equations  related  to  (i)  velocity,  (ii) 
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temperature  and  (iii)  the  remaining  equations.  For  example,  for  the  one  dimensional  planar  configuration, 
the  coordinates  are  chosen  such  that  the  wall  is  parallel  to  the  x  direction  and  y  is  the  direction 
perpendicular  to  the  wall.  The  velocity  in  the  x  direction  is  u  and  the  velocity  in  the  other  directions  is  zero. 
All  derivatives  in  the  x  direction,  except  pressure  gradient  in  Poiseuille  flow,  are  zero  everywhere  and 
mass  conservation  is  satisfied  automatically.  The  equations  involved  in  the  velocity  problem  from  the 
linearised  R26  moment  system  (LR26)  reduce  to  the  following  five  moment  system: 


du  dcrxy  _  Qp 
dt  dy  fix’ 


(70) 


dt 


+ 


d_\ 
dy  l 


_  2  _ 
u  +—qr 
5 


y 


(71) 


+  =  [?LL n 

dt  2  dy  3  V  2  A  **  dy  ’ 


(72) 


C  3  WL_  8  daxy 

ST  +  dy  2\  2  Amxy>’  5  dy  ’ 


(73) 


9RXy  |  dYxyy  __  7  L  -  14 

St  dy  6\  2  2  xy  5  dy 


(74) 


with  the  associated  wall  boundary  conditions  (59),  (63)  and  (65),  and  the  following  constitutive 
relationships: 


Qxyyy 


J!  xT/ 

7Zj  L  dy  ’  Vxyy 


_72_ JYz and  s  . /2 
35  Yx\  t:  L  dy  V  n  L  dy 


(75) 


6.0  NUMERICAL  SOLUTION  OL  THE  MOMENT  EQUATIONS 

The  moment  method  results  in  a  set  of  partial  differential  governing  equations  of  mixed  first  and  second 
order.  In  most  situations,  there  are  no  analytical  solutions  for  this  complex  set  of  governing  equations  and 
a  numerical  procedure  is  therefore  required.  However,  in  the  moment  equation  system,  the  momentum 
and  energy  conservation  equations  are  the  first  order  partial  differential  equations,  in  which  there  is  no 
explicit  gradient  transport  mechanism.  The  inadequacy  of  the  standard  finite  volume  method  for  the 
governing  equations  without  any  gradient  transport  terms  is  well  recognised  [22],  so  that  methods  for 
dealing  with  hyperbolic  problems  are  required.  In  the  case  of  low-speed  rarefied  gas  flow,  such  as  those 
found  in  micro-devices,  the  flow  is  parabolic  or  elliptic.  Using  a  hyperbolic  flow  solver  to  solve  elliptic 
flows  is  inefficient  and  expensive.  In  fact,  the  gradient  transport  mechanism  of  the  low  moments  is 
embedded  in  the  moments  one  order  higher.  For  example,  the  diffusion  of  u,  is  included  in  ay,  and  the 
diffusion  of  T  is  included  in  qu  as  reflected  by  the  underline  terms  on  the  right  hand  of  equations  (18)  and 
(19).  To  restore  the  gradient  transport  mechanism  back  to  the  moment  equation  system,  a  primitive 
variable  transformation  was  introduced  by  decomposing  the  moments  into  their  gradient  and  non-gradient 
components.  The  resultant  equations  of  the  non-gradient  components  have  the  general  convection- 
diffusion  form  [14,  15].  Alternatively,  the  gradient  transport  mechanism  can  be  introduced  into  the 
discritised  moment  equations  in  a  collocated  grid  by  a  properly  constructing  the  flux  in  the  control  volume 
face  [26]. 
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For  a  two  dimensional  steady  state  flow,  the  momentum  equation  (12)  can  be  written  in  Cartesian 
coordinates  (x,y)  as 


dpUU  |  dpVU  ^  o(Txx  |  da y  _  dp 

dx  dy  dx  dy  dx’ 

dpUV  dpVV  Saxy  do ^  _  dp 

dx  dy  dx  dy  dy  ’ 


(76) 

(77) 


in  which,  U  and  V  are  the  gas  velocity  inx  and  y  direction,  respectively.  In  a  collocated  grid  arrangement 
shown  in  Figure  3,  all  the  variables  are  stored  at  the  control  volume  centre  P  and  its  neighboring  points 
WW  W  P  E  EE  SS  S  N  NN.  A  discrete  approximation  to  equation  (76)  or  (77),  in  which  U  or  V  is  the 
primitive  variable,  is  achieved  by  integrating  the  equation  over  a  macro-control  volume  surrounding  the 
nodal  point  P, 


r  dpUU 
K  dx 


+ 


dx 


\ 

dxdy  + 

/» n 

C  e 

f dpVU  ,  SO 

dxdy  = 

c n 

Ce  f 

/ 

•J 

s  ^ 

w 

l  fy  dy  ) 

J 

s  ^ 

w  V 

dp 

dx 


dxdy 


(78) 


where  n,  s,  e  and  w  refer  to  the  location  of  the  space  average  of  any  quantity  prevailing  over  the  faces  of 
the  control  volume.  Using  the  mean-value  theorem,  equation  (78)  can  be  written  as: 


\pvu\-(puu\  K),-K). 


dx 


Sx 


dxdy  + 


(Pvu)-(Pru\ 


dy 


dy 


dxdy 


(79) 


6  dp ^ 

V  dx  jP 


dxdy. 


Figure  3:  Collocated  grid  arrangement. 

As  suggested  by  the  underline  terms  in  the  right  hand  side  of  equation  (18),  the  stress  in  the  momentum 
equation  is  proportional  to  the  velocity  gradient.  In  the  discritised  form  in  x  direction,  the  normal  stress  on 
the  nodal  points,  P,  W  and  E  should  have 


RTO-EN-AVT-194 


11  -15 


Application  of  the  Moment  Method  in  the 

Slip  and  Transition  Regime  for  Microfluidic  Flows 


ORGANIZATION 


K7 


cc  — 


A 


SU_ 

Sx 


(  ^ XX  )  II 


OC  — 


Jp 


A 


SU_ 

Sx 


and  (^)EX~ 


Jw 


su_ 

Sx 


(80) 


Je 


The  stress  at  the  control  volume  faces,  w  and  e  should  have  the  same  property  as: 

(O* 


oc  — 


(  su) 

i  /  \ 

(  su) 
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l  Sx)e 

(81) 


To  estimate  the  values  of  the  stresses  at  the  faces  of  the  control  volume,  it  is  important  to  preserve  this 
property.  To  guarantee  this,  it  is  convenient  to  interpolate  the  face  value  of  the  normal  stress  as: 
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Similar  interpolation  of  the  face  value  of  the  shear  stress  in  the  y  direction  can  be  obtained  as 
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Inserting  equations  (82)  and  (83)  into  equation  (79),  the  discritised  momentum  equation  in  the  x  direction 
can  be  expressed  by 


S(pUU)  5(pVU) 
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Equation  (84)  is  a  typical  discritised  convective  and  diffusion  equation  with  source  terms  on  the  right  hand 
side.  The  same  procedure  can  be  performed  to  equations  (77),  (13),  (18)  and  (19).  Numerical  methods  for 
solving  the  convective  and  diffusion  equations  are  well  documented  for  both  high  and  low  speed  flows 
[22],  For  the  convective  terms,  a  range  of  upwind  schemes  including  QUICK  [32],  SMART  [33],  and 
CUBISTA  [34]  are  available  in  the  literature.  The  diffusive  and  source  terms  are  discretised  by  a  central 
difference  scheme  and  the  CUBISTA  scheme  was  selected  for  the  present  study.  The  coupling  of  the 
velocity  and  pressure  fields  is  through  the  SIMPLE  [35]  or  PISO  [36]  algorithm.  For  steady  state  flow,  the 
system  of  equations  can  be  solved  iteratively  and  the  solution  procedure  is  summarised  as  follows: 

1)  Solve  Uj  at  iteration  n+ 1  using  the  values  of  other  variables  at  the  previous  iteration  n. 

2)  Solve  the  pressure  correction  equation  using  the  SIMPLE  or  PISO  algorithm  to  update  p  and  it,  at 
iteration  n+ 1. 
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3)  Solve  T,  a,j,  q„  mijk,  Ry,  A  at  iteration  n+ 1  using  updated  pressure  and  velocity  fields. 

4)  Update  the  boundary  conditions  according  to  equations  (36),  (37)  and  (39)-(46). 

5)  Return  to  step  1  and  repeat  until  residuals  of  each  governing  equation  reaches  a  specified 
convergence  criterion. 

Computationally,  it  is  more  expensive  to  solve  the  R26  equations  than  the  NSF  equations  but  this  cost  is 
necessary  to  capture  non-equilibrium  phenomena.  However,  the  advanced  computational  and  numerical 
techniques  developed  over  the  years  for  conventional  computational  fluid  dynamics  can  be  readily  adopted 
so  that  the  moment  method  can  be  applied  for  engineering  applications  where  non- equilibrium  effects  are 
important.  In  particular,  in  the  low-speed,  low  Reynolds  number  regime  where  it  can  be  costly  and 
difficult  to  get  meaningful  statistical  data  from  stochastic  methods,  the  R26  equations  have  a  significant 
computational  advantage. 

7.0  VALIDATIONS  AND  APPLICATIONS 

As  the  macroscopic  models  are  reduced  from  detailed  kinetic  theory,  some  kinetic  information  is  lost 
during  the  reduction  process.  Two  important  features  of  solutions  to  the  Boltzmann  equation  is  that  the 
distribution  function,  f  is  non- negative  and  the  solution  must  satisfy  the  //-theorem  [6].  However,  a 
solution  obtained  from  the  moment  equations  as  an  approximation  to / may  not  satisfy  these  constraints.  It 
is  proved  that  the  linearised  R13  equations  naturally  fulfil  the  //-theorem  [37].  For  full  R13  and  R26 
moment  equations,  numerical  analysis  is  required  to  assess  the  moment  realizability  and  the  //-theorem 
inequality  [14],  From  the  engineering  point  of  view,  it  is  necessary  to  validate  the  macroscopic  models 
against  available  kinetic  data  as  much  as  possible.  In  this  section,  we  will  first  use  the  linearised  moment 
equations  to  obtain  analytical  solutions  for  Kramers’  problem  and  Poiseuille  flow.  Numerical  techniques 
are  then  used  to  solve  the  moment  equations  for  complicated  flow  configurations.  The  analytical  and 
numerical  solutions  will  be  compared  with  kinetic  data  for  the  Boltzmann  equation  and  DSMC  data  to 
demonstrate  the  capability  as  well  as  the  limitation  of  the  macroscopic  modelling  approach. 

7.1  Kramers’  Problem  and  Knudsen  Layer 

The  Knudsen  layer  is  a  boundary  layer  that  bridges  the  wall  effect  and  the  bulk  flow.  Kramers’  problem 
[38],  first  formulated  in  1949,  is  the  most  basic  configuration,  where  the  effect  of  a  solid  wall  can  be 
investigated  without  additional  complications  found  in  more  realistic  geometries.  It  has  been  extensively 
studied  [39-43]  and  it  provides  a  useful  benchmark  for  the  development  of  macroscopic  models.  As  a 
half-space  problem,  it  can  be  readily  solved  by  the  linearised  equations  (70)-(75).  The  wall  is  set  at  y  =  0 

and  the  shear  stress,  a ,  is  constant  so  there  is  no  pressure  gradient  in  the  x  direction.  The  mean  free  path, 

A,  is  chosen  as  the  characteristic  length,  L,  so  that  A/L= 1.  Integration  of  equation  (71)  gives  the  velocity 
profile  as: 


(85) 


where  A  is  an  integration  constant  determined  by  the  wall  boundary  conditions.  The  supeiposition  of  the 
velocity  contributions  from  a  ,  qx  and  m  is  clearly  expressed  by  equation  (85).  This  is  an  important 

feature  that  is  not  present  in  the  NSF  equations,  in  which  qx  -  0  and  m  =  0  .  The  expressions  for  qx 

and  m  can  be  obtained  from  equations  (72)-(74)  as: 
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(86) 

(87) 


where  C \  and  C2  are  integration  constants.  In  the  above  solutions,  the  boundary  condition  that  as  y  — >  go  , 
qx,  mxyv  and  Rxy  will  remain  finite  has  been  used  to  remove  the  other  two  integration  constants  and  the 

terms  associated  with  them.  The  remaining  integration  constants,  A,  C\  and  C2,  are  determined  from  the 
wall  boundary  conditions  (59),  (63)  and  (65)  by 


c,=- 


C2  = 


(0.29  lcr  -t- 1 .658)  or 
0.101a2  +  1.128a  +  9.0 
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Inserting  equations  (86)-(90)  into  equation  (85),  the  final  expression  for  velocity  from  the  LR26  moment 
equations  reads 


y+2za(  Cie-{unfiJIy) 
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+  0. 178C2e 
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(91) 


In  the  linearised  R13  system  (LR13),  the  governing  differential  equations  (73)  and  (74)  are  replaced  by  the 
following  constitutive  relationships  [44]: 
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xyy 


16  2  [2 
15  L\Itt  8y 


and  R„, 


122  [2  dqx 
5  L\rc  ddy 


(92) 


along  with  the  boundary  conditions  (59)  and  (63)  without  the  higher  moments,  y/  and  Qx  .  For 
Kramers’  problem,  it  is  interesting  to  see  that  m  -  0  in  the  LR13  model.  The  velocity  field  for  the 
LR13  equations  is  readily  obtained  as 


 2-a 

(l3-2VHTr)«  +  4Vi(Tr 
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j  2  v 

y  + 
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(93) 


If  we  compare  equations  (91)  and  (93),  we  can  see  that  the  LR26  equations  provide  two  exponentials  to 
describe  the  Knudsen  layer  whilst  the  LR13  equations  only  contain  one  exponential. 

In  kinetic  theory,  the  defect  velocity  and  slip  coefficient  are  often  used  to  study  how  the  wall  affects  the 
velocity  profile  [39-43].  To  be  consistent  with  the  kinetic  solutions  obtained  from  the  Boltzmann 
equation,  a  reference  velocity  defined  by 
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is  used  to  scale  the  velocity  i.e. 
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A  defect  velocity,  Md,  is  defined  by  [40,  43] 


ud  -y  +  £  —  u 


(96) 


and  the  slip  coefficient,  is  determined  by 


hm  ud  =  0.  (97) 

y— >oo 

From  the  above  condition,  the  expression  of  £  for  the  LR26  moment  equations  can  be  written  as 
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and  the  defect  velocity  is  expressed  by 
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Similarly,  we  can  obtain  the  defect  velocity  and  slip  coefficient  for  the  LR13  equations,  respectively,  as 
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Figure  4  presents  the  analytical  solutions  of  the  defect  velocity  in  the  Knudsen  layer  from  the  moment 
equations  in  comparison  with  the  computational  results  from  the  three  kinetic  models  investigated  by 
Siewert  [43],  i.e.  the  BGK  model,  the  Williams  model  (the  collision  frequency  is  proportional  to  the 
magnitude  of  the  velocity),  and  the  hard  sphere  model.  For  the  case  of  a  —  0.9,  diffusive  reflection  from 
the  wall  dominates.  The  BGK  kinetic  model  produces  the  largest  defect  velocity  in  the  Knudsen  layer, 
particularly  close  to  the  wall  (y  <  0.5 ).  In  contrast,  the  results  from  the  Williams  and  hard  sphere  models 
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are  in  close  agreement  with  each  other.  The  analytical  solution  from  the  LR26  equations  generally  lies 
between  the  three  models  and,  beyond  y  =  0.5,  is  in  close  agreement  with  the  BGK  model.  However,  at 
the  wall,  all  kinetic  models  predict  a  higher  defect  velocity.  Conversely,  the  solution  obtained  from  the 
LR13  equations  undeipredicts  the  defect  velocity  significantly,  as  shown  in  Figure  4(a).  As  the  LR13 
system  involves  fewer  equations  and  boundary  conditions  than  the  LR26  system,  less  kinetic  information 
is  preserved  in  the  LR13  model.  Clearly  the  combination  of  two  exponentials  with  different  widths 
produces  an  improved  Knudsen  layer  velocity  profile.  It  is  expected  that  more  moments  and  their 
governing  equations  would  generate  more  exponentials  with  different  widths  to  recover  the  full  kinetic 
information.  As  the  value  of  the  accommodation  coefficient  decreases,  the  previous  observations  remain 
valid  although,  as  expected,  the  defect  velocity  increases,  as  shown  in  Figure  4(b).  In  this  case,  where 
a  =  0.1  ,  specular  wall  reflection  dominates.  The  results  shown  in  Figure  4  illustrate  that  the  LR26 
equations  work  well  for  walls  exhibiting  either  diffusive  or  specular  reflection. 


y  y 

Figure  4:  Defect  velocity  profile  for  Kramers’  problem:  comparison  between  the  moment  equation 
solutions  (lines)  and  kinetic  theory  (symbols)  [43].  (Note:  The  original  data  in  Ref  [43]  are 

presented  in  terms  of  the  mean  free  path  defined  by  l  =  {l/JxjA.. .  The  data  were 
converted  to  be  consistent  with  the  present  definition  of  the  mean  free  path,  A.). 

Equation  (85)  illustrates  that  the  velocity  in  the  Knudsen  layer  consists  of  contributions  from  a xy ,  qx ,  and 
m  .  For  the  LR13  equations,  mxyy  =0,  there  is  no  mechanism  for  m  to  contribute;  so  the  Knudsen 
layer  is  derived  solely  from  the  tangential  heat  flux,  qx.  In  the  NSF  system,  qx  and  mxyy  are  not  present, 

so  there  is  no  Knudsen  layer.  If  we  take  the  distance  to  the  wall  as  the  characteristic  length,  the 
corresponding  Knudsen  number  is  the  reciprocal  of  y~.  Figure  4  suggests  that  the  LR26  and  LR13 
equations  can  capture  the  Knudsen  layer  velocity  in  a  half  space  configuration  for  a  Knudsen  number 
equal  to  2  (y  -  0.5)  and  0.5  [y  =  2) ,  respectively.  However,  in  confined  geometries,  Knudsen  layers  from 

opposite  walls  will  overlap  [14]  and,  for  the  R26  and  R13  equations,  this  will  reduce  the  value  of  the 
Knudsen  number  to  1  and  0.25,  respectively.  These  results  provide  a  clear  indication  that  the  higher 
moment  equations  can  be  used  in  the  early  transition  regime  with  good  accuracy  to  account  for  the  wall 
effect  on  flows  in  MEMS. 
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7.2  Planar  Poiseuille  Flow  and  Knudsen  Minimum 

Channels  are  the  most  frequently  encountered  geometry  in  MEMS.  Planar  Poiseuille  flow  driven  by  a 
small  pressure  gradient  is  another  classic  case  often  used  to  check  the  validity  of  any  proposed 
macroscopic  models.  The  distance  between  the  plates  is  chosen  as  the  characteristic  length,  L,  so  that  the 
Knudsen  number, 


(102) 


The  two  plates  have  been  set  aty  =  ±L/2,  i.e.,  y  =  ±1/2 .  Symmetry  boundary  conditions  are  applied  at 
y  -  0,  and  the  flow  is  driven  by  a  constant  pressure  gradient 


dp 

dx 


=  -£• 


From  equation  (70),  the  shear  stress  can  be  expressed  by 

^  =  &  ■ 


(103) 


(104) 


Integration  of  equation  (71)  gives: 


u  =  - 


*  Z  -2  -  2 

—  y  ~m - qx+B 

2  2  Kn  5 


(105) 


in  which  B  is  an  integration  constant.  Again,  it  is  clearly  shown  from  equation  (105)  that  the  velocity 
profile  in  the  confined  geometry  has  contributions  from  not  only  axy,  but  also  from  mxyy  and  qx,  which  can 
be  obtained  from  equations  (72)-(74)  as: 
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where  C3  and  C4  are  the  integration  constants,  along  with  B,  are  determined  by  the  boundary  conditions 
(59),  (63)  and  (65)  as: 
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Substituting  equations  (106)-(108)  into  equation  (105),  the  velocity  profile  of  planar  Poiseuille  flow  is 
obtained  as 
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In  order  to  compare  with  the  linearized  Boltzmann  solution  [7,  45],  a  reference  velocity, 

u0=!;j2KT0 

is  used  to  renormalize  the  velocity  as 
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The  normalised  mass  flow  rate  of  planar  Poiseuille  flow,  Q,  can  be  obtained  by  the  integration  of  equation 
(113)  over  the  channel  width  as 
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Similarly,  we  can  obtain  the  velocity  profiles  from  the  LR13  equations  as 
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in  which  the  constant  C5  is  determined  from  the  boundary  conditions  as: 
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The  normalised  flow  rate  obtained  from  the  LR13  equations  is: 
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The  velocity  profiles  calculated  from  equations  (113)  and  (115)  from  the  LR26  and  LR13  moment 
equations,  respectively,  are  presented  in  Figure  5  at  a  range  of  Knudsen  numbers  in  comparison  with  the 
solution  from  the  linearised  Boltzmann  equation  [7].  At  Kn= 0.113,  which  is  just  beyond  the  slip-flow 
regime,  both  extended  hydrodynamic  models  predict  similar  values  of  velocity  and  are  all  close  to  the 
solution  obtained  from  the  Boltzmann  equation,  as  shown  in  Figure  5(a).  The  LR13  equations 
underpredict  while  the  R26  equations  slightly  overpredict  the  maximum  velocity.  Flowever,  as  the  value  of 
Kn  increases  and  the  flow  enters  the  transition  regime,  the  LR13  equations  overpredict  the  slip  velocity 
significantly.  In  contrast,  the  velocity  fields  predicted  by  the  LR26  equations  compare  very  well  to  the 
solution  obtained  from  the  Boltzmann  equation  for  both  Kn= 0.226  and  0.451,  as  shown  in  Figures  5(b) 
and  (c).  The  velocity-slip  predicted  by  the  LR26  equations  is  in  reasonable  agreement  with  the  value 
predicted  by  the  Boltzmann  equation  but  discrepancies  in  the  bulk  flow  begin  to  show  at  Kn= 0.903,  Figure 
5(e),  and  differ  at  A/i=1.128,  as  shown  in  Figure  5(f). 
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Figure  5:  Comparison  of  velocity  profiles  of  pressure-driven  Poiseuille  flow  at  different 
values  of  Knudsen  number  between  the  linearized  moment  equation 
solution  and  the  linearised  Boltzmann  equation  solution  [7], 

The  accurate  prediction  of  the  flow  rate  in  a  micro-channel  is  important  in  the  design  of  micro-devices.  To 
get  the  flow  rate  correct,  it  is  essential  that  the  predicted  velocity  profile  is  correct.  In  contrast,  the  correct 
prediction  of  the  flow  rate  cannot  guarantee  the  correct  velocity  profile.  Figure  6  shows  the  predicted  mass 
flow  rates  by  the  LR13  and  LR26  equations  in  comparison  to  the  solution  obtained  from  the  Boltzmann 
equation  [7,  44]  at  two  different  accommodation  coefficients.  When  Kn  <0.1  in  the  slip  regime,  the  flow 
rates  predicted  by  both  models  are  close  to  the  solution  obtained  from  the  Boltzmann  equation.  As  the 
value  of  Kn  increases,  non-equilibrium  effects  gradually  enter  the  central  flow  region,  the  LR26  equations 
follow  the  solution  obtained  from  the  Boltzmann  equation  reasonably  well  until  Kn  reaches  about  1.5.  As 
shown  in  Figure  6,  the  LR26  equations  predict  a  Knudsen  minimum  at  the  value  of  Kn  predicted  by  the 
Boltzmann  equation  for  both  a=  0.7  and  1.0.  The  R13  equations  are  close  to  the  Boltzmann  solution  up  to 
Kn  ~  0.4,  particularly  for  the  lower  value  of  a.  They  can  also  predict  a  Knudsen  minimum  but  at  a  value  of 
Kn  smaller  than  that  observed  by  the  Boltzmann  equation. 
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Kn 

Figure  6:  Comparison  of  mass  flow  rate  of  pressure-driven  Poiseuille  flow  at  different 
values  of  Knudsen  number  between  the  linearized  moment  equation 
solution  and  the  linearised  Boltzmann  equation  solution  [7,  44], 


One  of  the  interesting  effects  of  rarefaction  in  Poiseuille  flow  revealed  by  kinetic  theory  is  the  presence  of 
heat  flow  in  a  channel  without  any  temperature  gradient  [7].  This  phenomenon  is  also  captured  by  the 
LR26  moment  equations.  The  tangential  heat  flux  is  caused  by  the  pressure  gradient  and  the  Knudsen 
layer  contribution  as  expressed  by  the  first  and  second  terms  in  the  right-hand  side  of  equation  (107), 
respectively.  When  the  Knudsen  number  is  small,  heat  flows  in  the  opposite  direction  to  the  mass  flow  at 
the  centre  of  the  channel  but,  close  to  the  wall,  it  is  in  the  same  direction  as  the  mass  flow,  as  indicated  in 
Figure  7.  Flowever,  as  the  Knudsen  number  increases,  heat  flow  is  always  in  the  opposite  direction  to  the 
mass  flow.  This  is  a  high-order  rarefaction  effect  which  is  clearly  not  embedded  in  the  NSF  equations. 
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Figure  7:  Tangential  heat  flux  of  pressure-driven  Poiseuille  flow  at  different  values  of  Kn 
with  fully  diffusive  walls.  Lines:  the  LR26,  Symbols:  the  linearised  Boltzmann  equation 

[7].  qx  is  the  normalised  tangential  heat  flux  by  qx  =qx/{ P0 -J2RT  'j  =  qx/-Jl . 


7.3  Stokes’  Second  Problem  and  Oscillatory  Planar  Couette  Flow 

Many  MEMS  devices  contain  oscillating  parts  where  air  (viscous)  damping  plays  an  important  role.  To 
understand  the  damping  mechanisms,  it  is  essential  to  consider  non-equilibrium  or  rarefaction  effects.  In 
oscillatory  flows,  additional  non-equilibrium  effects  can  arise  if  there  are  an  insufficient  number  of 
intermolecular  collisions  during  one  cycle  of  the  oscillation.  To  quantify  the  extent  of  rarefaction  in  this 
respect,  Sharipov  and  Kalempa  [45,  46]  introduced  a  rarefaction  parameter,  9,  as 

0  =  -,  (118) 

co 

in  which  co  and  //  are  the  oscillation  and  the  intermolecular  collision  frequency,  respectively.  By  analogy 
with  the  spatial  Knudsen  number,  Kn,  which  is  based  upon  a  typical  length  scale,  it  is  convenient  to  adopt 
a  temporal  Knudsen  number,  Kn,,  which  is  defined  as  the  reciprocal  of  9,  to  express  the  extent  of  non¬ 
equilibrium  from  the  aspect  of  time  scale  [47],  i.e. 

Knt=^  =  ~.  (119) 

9  7J 

The  intermolecular  collision  frequency,  77,  can  be  estimated  by  pip  [6].  When  Kn,  «  1,  many 
intermolecular  collisions  occur  during  one  cycle  of  the  oscillation  so  that  unsteadiness  has  little  effect  on 
the  flow  to  approach  the  equilibrium  state.  On  the  other  hand,  in  the  regime  where  Kn,  »  1,  very  few 
intermolecular  collisions  occur  during  one  cycle  of  the  oscillation.  The  intermolecular  collisions  can  be 
neglected  and  the  free  molecular  formulation  of  the  Boltzmann  equation  can  be  used.  In  a  system  where 
Kn,  ~  1,  the  number  of  the  intermolecular  collisions  is  not  sufficiently  large  for  the  flow  to  reach  an 
equilibrium  state  during  one  oscillation  cycle.  Rarefaction  effects  have  to  be  taken  into  account  when 
modeling  flows  in  this  regime. 
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Analytical  solutions  of  the  Navier-Stokes  equation  for  Stokes’  second  problem  and  oscillatory  planar 
Couette  flow,  with  either  first-order  or  second-order  velocity-slip  boundary  conditions,  can  be  found  in  the 
works  of  Sharipov  and  Kalempa  [45,  46],  Park  et  al.  [48]  and  Hadjiconstantinou  [49].  Emerson  et  al.  [50] 
studied  nonplanar  oscillatory  Couette  flow  in  the  entire  range  of  Knudsen  number.  The  analytical  solutions 
for  the  LR13  and  LR26  moment  equations  can  also  be  obtained;  they  are  a  linear  combination  of  4  and  6 
exponentials  with  different  widths  for  the  LR13  and  LR26  equations,  respectively.  However,  the 
expressions  for  the  coefficients  and  the  widths  are  rather  complicated  and  cumbersome.  It  is  therefore 
more  convenient  to  solve  the  one  dimensional  LR13  and  LR26  equations  numerically  [47].  Equations  (70) 
-(74)  can  readily  be  solved  using  a  central  difference  scheme  to  discretize  the  spatial  derivatives  and  a 
second-order  Crank-Nicholson  method  for  the  temporal  terms. 

Stokes’  second  problem  is  a  half  space  configuration  with  the  oscillating  wall  located  at  y  =  0.  The  mean 
free  path,  A,  is  chosen  as  the  characteristic  length  so  that  A/L  in  equations  (70)  -  (75)  is  equal  to  1.  The  gas 
molecules  are  assumed  to  diffusively  reflect  from  the  oscillating  wall  and  the  gradients  of  all  variables  are 
equal  to  zero  at  the  open  end  of  the  solution  domain.  To  solve  the  system  of  equations,  200  to  400  equi- 
spaced  grid  points  are  used  depending  on  the  oscillation  frequency.  For  oscillatory  planar  Couette  flow, 
the  lower  plate  is  located  at  y  =  0  and  remains  stationary.  The  upper  plate,  (y  =  L),  oscillates  at  a  fixed 
frequency.  The  gap  between  the  two  plates,  L,  is  chosen  as  the  characteristic  length  scale,  so  that  A/L  in 
equations  (70)  -  (75)  is  equal  to  Kn.  The  flow  is  solved  with  200  grid  points  across  the  one-dimensional 
domain.  For  both  configurations,  the  moving  wall  oscillates  sinusoidally  according  to 

ua  =  uw  sin(®0,  (120) 

where  a  and  uw  are  the  oscillating  frequency  and  magnitude  of  the  wall,  respectively,  and  u0  is  the 
instantaneous  wall  velocity.  The  numerical  time  interval  is  taken  as  one  hundredth  of  a  period.  The 
accommodation  coefficients  of  the  lower  and  upper  walls  are  assigned  a  value  of  unity.  All  computations 
are  impulsively  started  with  zero  velocity  everywhere  and  50  to  100  periods  are  computed  to  remove  the 
quiescent  initial  effects  and  to  ensure  that  a  quasi-steady  periodic  state  is  reached.  Computed  results  for 
the  macroscopic  models  are  then  compared  with  data  from  kinetic  theory  to  assess  their  validity. 

7.3.1  Stokes’  Second  Problem 

The  dynamic  velocity  profiles  from  the  linearised  moment  equations  at  four  points  in  time,  corresponding 
to  cot  =  0,  n/ 2,  n,  In/2,  are  presented  in  Figure  8.  At  the  lower  value  of  Kn,  equal  to  0.1,  where  many 
intermolecular  collisions  occur  during  one  cycle  of  oscillation,  the  results  from  both  extended  continuum 
models  are  nearly  identical,  as  shown  in  Figure  8(a).  Analogous  to  the  classification  for  steady  flow,  we 
can  denote  the  regime  of  Kn,  <0.1  as  the  hydrodynamic  regime.  As  Kn,  increases  to  0.5,  in  Figure  8(b), 
the  solutions  from  the  LR13  and  LR26  equations  start  to  deviate  from  each  other.  At  Kn,  equal  to  unity, 
well  into  the  transition  regime,  the  velocity  profiles  from  the  LR13  equations  are  significantly  different 
from  those  of  the  LR26  equations  not  only  in  the  region  close  to  the  wall  but  also  in  the  bulk  flow,  as 
indicated  in  Figure  8(c).  The  LR13  equations  predict  a  larger  slip  velocity  than  the  LR26  equations  which 
is  due  to  the  LR13  equations  lacking  a  proper  mechanism  to  capture  the  Knudsen  layer  velocity  profile 
near  the  wall  [2 1  ]  and  the  oscillation  aggravates  the  situation.  In  general,  the  slip  velocity  increases  as  Knt 
increases  so  that  the  velocity  on  the  wall  reduces. 


RTO-EN-AVT-194 


11  -27 


Application  of  the  Moment  Method  in  the 

Slip  and  Transition  Regime  for  Microfluidic  Flows 


ORGANIZATION 


0.6 

0.4 

0.2 

0.0 


-0.2 


-0.4 


-0.6 


y 


Figure  8:  Dynamic  velocity  profiles  of  Stokes’  second  problem  from  the  extended 
hydrodynamic  models  at  four  different  times.  Adapted  from  Ref.  [47]. 


The  predicted  values  of  the  velocity  amplitude  at  the  wall  from  each  of  the  hydrodynamic  models  are 
plotted  against  Kn ,  in  Figure  9(a)  and  compared  to  the  kinetic  data  from  the  linearised  BGK  model  [45]. 
The  Navier-Stokes  equations  always  under-predict  the  velocity  amplitude  on  the  wall  even  in  the 
hydrodynamic  region.  Both  extended  macroscopic  models  improve  the  prediction  of  um  on  the  wall 
significantly,  as  shown  in  Figure  9(a),  up  to  Knt  =  0.5  with  the  LR26  equations  showing  better  agreement 
than  the  LR13  equations. 


Figure  9:  Predicted  velocity  and  shear  stress  amplitudes  on  the  wall  of  Stokes’ 
second  problem  from  macroscopic  models  in  comparison  with 
the  data  from  the  linearised  BGK  kinetic  equation  [45]. 


To  compare  with  kinetic  data  [45],  the  shear  stress,  axy,  is  renormalised  as 

PRTq 

2Po  Uw 


(121) 


The  predicted  value  of  the  shear  stress  amplitude  on  the  wall,  IT„„  from  the  hydrodynamic  models  and 
kinetic  theory  are  presented  in  Figure  9(b).  In  the  hydrodynamic  regime,  where  Kn,  <  0.05,  all  the 
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hydrodynamic  models  predict  similar  values  of  1T„,  on  the  wall  and  are  in  good  agreement  with  kinetic 
theory.  When  Hm  is  exceeds  0.1,  the  NS  equations  start  to  underestimate  1 1„„  while  both  the  LR13  and 
LR26  equations  are  still  able  to  follow  the  kinetic  data.  The  LR13  equations  tend  to  oveipredict  Wm  when 
Kn,  >  0.5. 

The  corresponding  velocity  and  shear  stress  phases  on  the  wall,  <pu  and  <pr,  as  defined  by  Sharipov  and 
Kalempa  [45],  are  presented  in  Figure  10.  The  velocity  phase,  (pu ,  has  a  maximum  value  in  the  transition 
regime  as  shown  in  Figure  10(a).  The  NS  and  LR13  equations  fail  to  predict  any  maximum  value.  In 
contrast,  the  LR26  equations  predict  a  maximum  value  of  (pu ,  but  the  value  is  larger  than  that  from  kinetic 
theory.  In  general,  the  agreements  of  <pu  between  the  hydrodynamic  models  and  kinetic  theory  are  not 
satisfactory  in  the  slip  and  transition  regimes,  although  the  LR26  model  improves  the  prediction 
substantially.  Flowever,  it  should  be  pointed  out  that  the  maximum  value  of  <pu  on  the  wall  is  less  than  4% 
of  one  period.  Conversely,  all  extended  hydrodynamic  models  predict  the  shear  stress  phase,  <pp ,  very  well 
up  to  Kn,  =  1,  as  illustrated  in  Figure  10(b).  The  NS  equations  underestimate  the  value  of  tpp,  from  a  very 
low  value  of  Kn,. 


Figure  10:  Predicted  velocity  and  shear  stress  phases  on  the  wall  of  Stokes’ 
second  problem  from  macroscopic  models  in  comparison  with  the  data 
from  the  linearised  BGK  kinetic  equation  [45],  Adapted  from  Ref.  [47], 


7.3.2  Oscillatory  Planar  Couette  Flow 

For  oscillatory  planar  Couette  flow  under  rarefied  conditions,  non-equilibrium  effects  will  arise  from  both 
the  separation  distance  and  the  oscillation  of  the  walls.  Therefore,  both  Knudsen  numbers,  Kn  and  Knt, 
defined  upon  a  length  and  time  scale,  respectively,  are  required  to  measure  the  extent  of  non-equilibrium. 
Traditionally,  the  Stokes  number,  S,  defined  by 


5  = 


(122) 


is  often  used  to  characterize  the  balance  between  the  unsteady  and  viscous  effects  [48,  49].  To  reflect  non¬ 
equilibrium  effects  due  to  the  oscillation,  the  temporal  Knudsen  number,  Kn,,  needs  to  be  assessed  to 
measure  the  oscillatory  effect  on  the  flow,  which  is  related  to  Kn  and  S  through 
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Knt  =  —(SKn )~ . 
n 


(123) 


Equation  (123)  indicates  that  Kn,  is  proportional  to  the  square  of  the  product  of  Kn  and  S.  Increasing  either 
Kn  or  S  can  lead  to  a  significant  increase  in  the  value  of  Kn,. 

Figure  1 1  shows  the  velocity  amplitude,  um,  between  the  two  walls  for  oscillatory  planar  Couette  flow. 
Results  at  different  values  of  Kn  and  Kn,  in  the  early  transition  regime  are  compared  with  the  DSMC  data 
[48]  and  the  solution  of  the  linearised  BGK  equation  [46].  For  Kn  <  0.1,  the  velocity  amplitudes  from  the 
LR13  and  LR26  equations  agree  with  both  DSMC  data  and  kinetic  theory  up  to  Kn,  =1.  When  Kn  =  0.5 
and  Kn,  =  0.16,  the  LR26  equations  are  in  good  agreement  with  the  DSMC  data.  Flowever,  when  Kn  =  1 
and  Kn,  =  0.64,  the  LR26  equations  can  capture  the  velocity  amplitude  reasonably  well  but  the  LR13 
equations  cannot  follow  the  DSMC  data  in  the  whole  domain,  as  illustrated  in  Figure  1 1(a).  For  the  case 
with  Kn  =  0.886  and  Kn,  =1,  both  the  LR26  and  LR13  equations  fail  to  capture  the  velocity  amplitude 
correctly,  as  indicated  in  Figure  1 1(b). 
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Figure  11:  Normalised  velocity  amplitude  of  oscillatory  planar  Couette  flow.  Lines  are  from  the 
macroscopic  models,  (a)  Symbols  are  the  DSMC  data  are  digitised  from  Park  et  al.  [48].  (b)  Symbols 
are  the  linearised  BGK  data  are  digitised  from  Sharipov  and  Kalempa  [46].  Adapted  from  Ref.  [47]. 


One  of  the  unique  features  in  one-dimensional  steady-state  Couette  flow  is  the  constant  shear  stress  across 
the  whole  domain  at  any  value  of  Kn,  even  inside  the  Knudsen  layer.  However,  this  feature  no  longer 
exists  when  one  of  the  plates  oscillates.  The  amplitude  of  the  shear  stress  is  therefore  presented  in  Figure 
12.  When  both  Kn  and  Kn,  are  small,  the  extended  hydrodynamic  models  agree  with  kinetic  theory,  but 
when  these  Knudsen  numbers  are  large,  the  discrepancies  between  them  are  very  noticeable.  The  complex 
interplay  of  the  length  and  time  scales  exacerbates  the  capabilities  of  the  macroscopic  equations  to  capture 
the  non-equilibrium  effects  at  high  values  of  either  Kn  or  Kn,. 
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Figure  12:  Normalised  shear  stress  amplitude  of  the  oscillatory  planar  Couette  flow. 
Lines  are  from  the  macroscopic  models.  Symbols  are  the  linearised  BGK  data 
are  digitised  from  Sharipov  and  Kalempa  [46].  Adapted  from  Ref.  [47]. 


The  velocity  and  shear  stress  phases,  (pu  and  (pP,  of  oscillatory  planar  Couette  flow  are  presented  in  Figure 
13.  When  Kn  =  0.0886  and  Kn,  =  0.1,  the  Stokes  number,  S,  is  equal  to  4.473.  The  LR13  and  LR26 
equations  predict  both  cpu  and  cpP  in  excellent  agreement  with  kinetic  theory  [46].  For  Kn  =  0.0886  and  Kn, 
=  1,  the  Stokes  number  reduces  to  1.415.  Although  the  LR26  equations  are  not  able  to  predict  um  and  nm 
very  well  for  this  particular  case,  as  shown  in  Figures  1 1(b)  and  12,  they  do  capture  the  phases,  <pu  and  <pP, 
quite  accurately  in  comparison  with  the  kinetic  data,  as  indicated  in  Figure  13.  In  contrast,  when  Kn  = 
0.0886  and  Kn,  =  1,  both  the  LR13  and  LR26  equations  are  able  to  predict  the  velocity  and  shear  stress 
amplitudes  correctly,  as  shown  in  Figure  11(b),  but  fail  to  predict  the  corresponding  phases.  The  LR26 
equations  manage  to  follow  the  kinetic  data  within  one  quarter  of  L  away  from  the  oscillating  wall  but 
then  start  to  overpredict  (pu  and  cpP  as  we  move  further  away  from  the  moving  wall,  as  indicated  in  Figure 
13.  The  Stokes  number  for  this  case  is  very  high,  equal  to  14.14.  At  such  a  large  Stokes  number,  inertia 
plays  a  major  role  in  the  flow  dynamics.  An  analysis  of  the  velocity  amplitudes  shown  in  Figure  11(b) 
indicates  that  the  penetration  depth  for  this  case  (Kn  =  0.0886,  Kn,  =  1)  is  0.38T  away  from  the  oscillating 
wall.  The  lower  stationary  wall  hardly  feels  the  any  effect  from  the  motion  of  the  upper  wall.  It  is 
interesting  to  see  that  the  LR26  equations  can  predict  <p„  and  (pP,  within  the  Stokes  layer  while  the  LR13 
equations  fail  to  follow  the  trend  of  (pu  and  <pP,  throughout  the  whole  domain. 
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Figure  13:  Velocity  and  shear  stress  phases  of  oscillatory  planar  Couette  flow.  Lines  are 
from  the  macroscopic  models.  Symbols  are  the  linearised  BGK  data  digitised 
from  Sharipov  and  Kalempa  [46].  Adapted  from  Ref.  [47]. 


7.4  Micro  Driven  Cavity  Flow 

The  driven  cavity  problem  is  a  well  known  benchmark  problem  for  testing  and  verifying  continuum 
solvers  [51].  The  problem  geometry  is  simple  and  two-dimensional,  yet  the  flow  pattern  and  heat  transfer 
in  the  cavity  are  encountered  in  many  engineering  applications.  However,  few  studies  [52-54]  have  been 
carried  out  to  examine  the  mass  and  heat  flow  patterns  in  a  driven  cavity  for  non-equilibrium  gas  flows. 
Recently,  John  et  al.  [55]  performed  a  series  of  DSMC  simulation  of  gas  flow  in  a  micro  lid-driven  cavity 
for  a  range  of  Knudsen  number,  Reynolds  number  and  Mach  number.  The  simulation  results  reveal 
interesting  non-equilibrium  phenomena  in  the  cavity  and  provide  useful  data  to  benchmark  the  extended 
hydrodynamic  equations. 

The  configuration  of  the  square  driven  cavity,  of  size  L,  is  shown  in  Figure  14.  The  notations  A,  B,  C,  and 
D  shown  in  the  figure  denote  the  four  comers  of  the  cavity.  The  top  lid  moves  with  a  fixed  tangential 
velocity  Uw  in  the  positive  x  direction,  and  the  other  walls  are  stationary.  The  wall  temperature  is  set  to  the 
reference  temperature,  i.e.,  Tw  =  273  K.  Any  variation  in  Knudsen  number  is  achieved  by  changing  the 
density  conditions  in  the  cavity,  i.e.,  the  reference  pressure pa.  The  R26  moment  equations  discritised  on  a 
129x129  uniform  grid  are  numerically  solved  for  a  range  of  Knudsen  numbers  in  the  slip  and  early 
transition  regimes  in  comparison  with  the  DSMC  data. 
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Figure  14:  Configuration  of  lid-driven  cavity. 


One  of  the  interesting  phenomena  in  a  driven  cavity  under  non-equilibrium  condition  is  the  gas 
temperature  distribution  and  heat  flux  direction.  Shown  in  Figure  15(a)  is  the  DSMC  simulation  of  heat 
flux  stream  traces  overlaid  on  temperature  contours  at  Kn  =  0.2  with  Uw  =10  m/s.  From  the  temperature 
contours,  a  cold  region  is  found  toward  the  left  comer  of  the  cavity,  whereas  a  hot  region  is  observed 
toward  the  right  comer  of  the  cavity.  More  interestingly,  it  is  noted  that  the  direction  of  heat  flow  is 
generally  from  the  cold  to  the  hot  region,  as  illustrated  by  the  heat  flux  streamlines.  This  represents  a 
counter-gradient  heat  flux  made  possible  by  the  rarefied  flow  conditions.  Under  non-equilibrium  flow 
conditions,  various  factors  such  as  expansion  cooling,  viscous  heat  generation,  compressibility,  and 
thermal  creep  could  significantly  affect  flow  and  heat  transfer.  For  the  driven  cavity  case,  an  expansion 
cooling  (gas  temperature  T  less  than  wall  temperature  Tw)  occurs  at  the  top  left  comer  of  the  cavity  due  to 
a  drop  in  pressure  which  results  in  heat  transfer  from  the  wall  to  the  gas,  whereas  viscous  heat  generation 
(r  >  Tw)  results  in  heat  transfer  from  the  gas  to  the  wall  toward  the  right  corner  of  the  cavity.  The  direction 
of  heat  transfer  is  governed  by  both  expansion  cooling  and  viscous  dissipation.  The  heat  transfer  plots  in 
the  driven  cavity  indicate  that  thermal  energy  transfer  need  not  always  be  from  a  hot  region  to  a  cold 
region  as  continuum  theory  dictates.  The  temperature  field  and  the  heat  flux  streams  predicted  by  the  R26 
equations  are  presented  in  Figure  15(b).  In  comparison  with  Figure  15(a),  the  R26  equations  capture  the 
temperature  field  and  the  heat  flux  streams  quite  well  except  that  the  R26  equations  overpredict  the 
temperature  drop  and  rise  in  the  top  left  and  right  comer  of  the  cavity,  respectively.  On  the  other  hand,  the 
NSF  equations  cannot  predict  thermal  field  neither  quantitatively  nor  qualitatively,  as  illustrated  in  Figure 
15(c). 
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Figure  15:  Comparison  of  heat  flux  stream  traces  overlaid  on  temperature  contours  computed 
by  (a)  DSMC,  (b)  R26  and  (c)  NSF  for  the  driven  cavity  at  Kn  =  0.2  and  Uw  =  10  m/s. 


The  predicted  profiles  of  velocity,  u  and  v,  shear  stress,  axy,  and  temperature,  T,  by  the  R26  moment 
equations,  along  the  vertical  and  horizontal  lines  across  the  cavity  centre  are  plotted  in  Figures  16-18  in 
comparison  with  the  DSMC  data  [55]  of  Uw  =  50  m/s  and  Kn  =  0.1  and  0.5.  When  Kn  =  0.1,  the  predicted 
velocities  in  both  directions  are  in  good  agreement  with  the  DSMC  data  as  shown  in  Figure  16.  As  Kn 
increases  to  0.5,  the  DSMC  results  show  that  the  values  of  the  slip  velocity  on  the  walls  increase.  Flowever, 
the  R26  equations  overpredict  about  1 0%  of  the  slip  velocity  on  the  top  moving  wall.  Away  from  the  lid, 
the  predicted  u  velocity  agrees  with  the  DSMC  data  quite  well  as  shown  in  Figure  16(a).  The  agreement 
with  the  DSMC  data  for  v  velocity  at  Kn  =  0.5  is  not  good  close  to  the  walls  as  well  as  away  from  the 
walls  as  plotted  in  Figure  16(b).  Flowever,  the  R26  equations  predict  the  values  of  shear  stress  quite 
accurately  on  the  vertical  line  for  both  Knudsen  numbers  as  indicated  in  Figure  17(a).  Figure  17(b)  shows 
that  the  value  of  shear  stress  is  overpredicted  by  the  R26  equations  on  the  horizontal  line  for  Kn  =  0.5. 
Flere  a„  =  /AJJL  is  a  reference  stress  for  normalization. 


u/Uw  x  /L 

Figure  16:  Comparison  of  normalised  velocity  profiles  across  the  cavity  centre. 
Lid  velocity  Uw=  50  m/s.  Symbols:  DSMC  [55];  Lines:  the  R26  equations. 

(a)  u  velocity  on  a  vertical  line;  (b)  v  velocity  on  a  horizontal  line. 
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Figure  17:  Comparison  of  normalised  shear  stress  profiles  across  the  cavity  centre.  Lid  velocity 
Uw=  50  m/s  and  cr0  =  //Uw/Z_.  Symbols:  DSMC  [55];  Lines:  the  R26  equations. 

(a)  Shear  stress  on  a  vertical  line;  (b)  Shear  stress  on  a  horizontal  line. 


The  predicted  temperature  profiles  on  the  vertical  and  horizontal  lines  across  the  cavity  are  plotted  in 
Figure  18  in  comparison  with  the  DSMC  data,  which  are  scattered.  It  is  computationally  expensive  to 
reduce  statistical  noise  for  temperature  field  with  DSMC  simulation,  particularly  for  flows  in  MEMS. 
Figure  18(a)  shows  that  the  gas  temperature  close  to  the  moving  lid  is  higher  than  the  lid  temperature.  At 
Kn  =  0.1,  the  temperature  of  gas  increases  as  it  moves  into  the  cavity,  then  drops  after  it  enters  into  more 
than  10%  of  the  cavity  depth.  The  R26  equations  capture  this  trend  well  in  good  agreement  with  the 
DSMC  data.  At  Kn  =  0.5,  the  maximum  temperature  of  the  gas  on  the  vertical  line  across  the  cavity  centre 
is  on  the  moving  wall.  The  R26  equations  underpredict  the  temperature  jump  significantly.  The 
temperature  of  the  gas  is  lower  than  Tw  on  the  left  while  higher  than  Tw  on  the  right  as  shown  in  Figure 
18(b).  The  predictions  of  the  R26  equations  are  in  reasonable  agreement  with  the  DSMC  data. 


Figure  18:  Comparison  of  normalised  temperature  profiles  across  the  cavity  centre. 
Lid  velocity  Uw=  50  m/s.  Symbols:  DSMC  [55];  Lines:  the  R26  equations. 

(a)  temperature  on  a  vertical  line;  (b)  temperature  on  a  horizontal  line. 
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7.5  Thermal  Transpiration  and  Knudsen  Pump 

The  phenomenon  of  thermal  transpiration  or  creep  was  discovered  by  Reynolds  in  1879,  in  which  a  gas 
will  move  along  a  solid  surface  due  to  inequalities  of  temperature  [56].  It  has  been  investigated 
theoretically  and  experimentally  by  various  researchers  over  more  than  a  century  [57-65].  It  is  assumed 
that  equal  numbers  of  molecules  arrive  at  the  wall  from  the  hot  and  cold  regions.  Molecules  arriving  from 
the  hot  region  will  have,  on  average,  a  higher  velocity  than  those  arriving  from  the  cold  region.  Since  the 
molecules  are  reflected  diffusively  at  the  wall,  the  resultant  force  on  the  wall  due  to  the  molecular 
collisions  acts  towards  the  cold  region.  An  equal  and  opposite  force  is  felt  by  the  gas  molecules,  giving 
rise  to  a  flow  towards  the  hot  region.  Once  the  fluid  starts  to  creep  along  the  wall,  the  moving  fluid  layer 
interacts  with  the  stagnant  fluid  layers  adjacent  to  it,  inducing  a  boundary  layer.  This  phenomenon  was 
used  by  Knudsen  in  1910  to  construct  the  first  multistage  thermal  transpiration  pump  [58,  59]. 

A  microscale  gas  pump  is  often  required  to  form  a  complete  micro  system.  Many  issues  have  been 
encountered  in  attempting  to  shrink  full-scale  conventional  pumps  to  microscales,  such  as  manufacturing 
tolerances,  pump  oil,  thermal  inefficiency  and  short  life  time.  The  Knudsen  pump  has  the  advantages  of  no 
moving  parts  and  supplementary  pumping  fluids.  A  typical  multistage  Knudsen  pump  is  a  long  pipe  or 
channel  with  a  periodic  structure  consisting  of  alternately  arranged  narrow  and  wide  dimensions.  The 
temperature  along  the  pipe  or  channel  is  also  periodic  with  a  distribution  increasing  in  the  narrow  parts 
and  decreasing  in  the  wide  parts.  In  the  narrow  parts  where  temperature  increases,  the  pressure  increases 
due  to  the  non-equilibrium  phenomenon  of  thermal  creep.  In  the  wide  section,  the  thermal  creep  effect  is 
less  profound  and  the  gas  flow  is  closer  to  the  continuum  regime.  One  of  the  modem  versions  of  the 
Knudsen  pump  was  designed  by  Pham-Van-Diep  et  al.  [66],  the  7th  stage  of  which  is  illustrated  in  Figure 
19.  The  performance  of  the  modem  Pham-Van-Diep  pump  has  been  studied  by  Vargo  et  al.  [67]  and 
Muntz  et  al.  [68]  using  kinetic  theory.  Other  modem  Knudsen  pumps  are  reported  by  McNamara  and 
Gianchandani  [69]  and  Pharas  and  McNamara  [70]. 
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Figure  19:  Illustrative  ith  stage  of  a  Knudsen  pump.  Adapted  from  Ref.  [67], 


One  of  the  major  challenges  in  the  design  of  the  Knudsen  pump  is  to  ensure  good  thermal  isolation 
between  the  hot  and  cold  ends  in  the  alternating  heating  and  cooling  structure.  One  way  to  avoid  this 
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difficulty  is  to  replace  the  alternating  heating  system  by  a  central  heating  system.  Illustrated  in  Figure  20 
is  a  3  staged  Knudsen  pump  heated  by  a  central  heater  on  the  right  side  to  generate  the  necessary 
temperature  gradient  for  thermal  creep.  The  temperature  gradient  is  largely  determined  by  the  thermal 
conductivity  of  the  channel  material  and  the  length  of  the  channels.  As  an  example,  it  is  assumed  that 
temperature  is  linearly  distributed  between  the  hot  and  cold  end  with  a  temperature  difference  AT.  The 
width  of  the  narrow  channel  is  designed  to  be  1  pm.  The  width  of  the  wide  channel  is  5  times  that  of  the 
narrow  ones.  The  R26  equations  are  solved  for  different  values  of  AT  and  Knudsen  number  Kna,  which  is 
based  on  the  inlet  condition. 


Figure  20:  An  illustration  of  a  central  3  stage  heated  Knudsen  pump  and  gas  flow  velocity  vectors  in  the 
pump.  Red  colour  indicates  hot  and  blue  cold.  The  temperature  on  the  walls  is  linearly  distributed. 


Shown  in  Figure  20  are  the  velocity  vectors  in  the  central  heated  Knudsen  pump  for  the  case  of  Kna  =  0.1 
and  AT  =  2 K.  Cold  gas  is  driven  towards  the  hot  end  in  the  narrow  channels  by  thermal  transpiration.  In 
the  wide  channels,  only  the  gas  close  to  the  walls  is  creeping  towards  the  hot  end.  The  gas  away  from  the 
walls  in  the  wide  channels  is  pushed  towards  the  cold  end  by  the  pressure  generated  in  the  narrow 
channels.  The  pressure  changes  Ap  along  the  channel  centre  line  from  the  inlet  to  the  outlet  are  plotted  in 
Figure  21  for  Kna  =0.1  and  AT  =  IK,  2 K  and  3 K,  respectively.  The  pressure  is  built  up  in  the  narrow 
channels  and  drops  in  the  wide  channels  and  when  the  flow  changes  direction.  The  width  ratio  of  the 
narrow  and  wide  channels  will  affect  the  pressure  distribution.  When  the  temperature  difference  between 
the  hot  and  cold  end  increases,  the  pressure  difference  between  the  inlet  and  outlet  increases. 
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Figure  21 :  Pressure  change  from  the  inlet  to  the  outlet  along  the  channel  centre  line. 


8.0  SUMMARY 

Gas  flows  in  micro-electro-mechanical-systems  suffer  from  non-equilibrium  effects.  The  conventional 
hydrodynamic  model,  the  Navier-Stokes-Fourier  equations,  is  unable  to  capture  these  effects  correctly  in 
the  slip  and  early  transition  regimes.  In  this  lecture,  we  demonstrate  that  hydrodynamic  models  can  be 
extended  to  account  for  non-equilibrium  effects  by  introducing  high  moment  equations,  which  can  be 
derived  from  kinetic  theory,  into  hydrodynamic  models.  From  Kramers’  problem  to  the  flow  in  a  lid 
driven  cavity,  it  is  indicated  that  the  R26  moment  equation  model,  based  on  Grad’s  moment  method  and 
Struchtrup  and  Torrilhon’s  regularization  procedure,  can  be  used  to  predict  non-equilibrium  gas  flows 
fairly  accurately  up  to  a  Knudsen  number  of  0.5.  The  conventional  numerical  techniques  for  low  speed 
flow  in  confined  geometries  can  be  readily  used  to  solve  the  moment  equations  for  engineering 
applications  in  MEMS. 
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